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Abstract 

The theory of the resonance radiation force is studied 
as it applies to the slowing of a two-level atom using a 
swept frequency laser. Ihe theory is developed in detail 
for this case. A single partial differential equation called 
the Fokker-Planck equation is found to describe the atomic 
motion for many cases. This equation is solved numerically 
for a one dimensional geometry. Results from this numerical 
solution are found to be within 10?J of experimental results. 
Analysis of this problem in one dimension leads to the 
prediction that for fixed amounts of frequency sweeping, a 
faster scan rate will yield smaller full width at half max¬ 
imum (FVfHIv!) spread in the velocity distribution for the atomic 
beam. This program also predicts that the atoms can be slowed 
arbitrarily close to 0 m/sec using this technique. The 
program used for this analysis is included. 





















LASER COOLING OF NEUTRAL ATOMS 


I. Introduction 


Background 

The theory of interaction of light with neutral atoms 
was first developed by Einstein (Ref 1) in 1917. In this paper 
Einstein introduced the concept of coefficients that describe 
absorbtion, spontaneous and stimulated emmission. He described 
how these coefficients allow derivation of the Planck radiation 
law and developed the theory of motion of atoms when interact¬ 
ing with light. The momentum transferred to the atoms, and the 
fluctuations in this transferred momentum, were calculated 
using quantized energy theory and the idea of the canonical 
distribution of states in thermal radiation. The first exper¬ 
imental work showing transfer of momentum to an atomic beam 
was performed by Frisch (Ref 2). Using conventional light 
sources, Frisch deflected an atomic beam. The advent of the 
laser, with high intensities and monochromatic nature, greatly 
enhanced the experimental worlds ability to slow atoms using 
light. 

In 1975 Hansch and Schawlow (Ref 3) first proposed that 
a standing wave laser could slow atoms using the spontaneous 
force. Very slow *« / 0 atomic beams are useful for 

atomic clocks (Ref 4 and 5)» doppler free spectroscopy (Ref 6), 
and insertion into single atom light traps (Ref 7). 





























In 1979 Balykin, Letokhov, and Minogin reported exper¬ 
imentally slowing atoms using a traveling wave constant fre¬ 


quency laser to slow atoms. The reduction in velocity was 
small* but the distribution functions full width at half max¬ 
imum (FWHM) was substantially reduced (Ref 8). In 1982 Phillips 
and Metcalf reported the first observation of slowing using a 
fixed frequency laser. The difference between the two methods 
is how the laser stays in resonance with the atomic beam. In 
the first case the frequency is changed to stay in resonance, 
and in the second a magnetic field is used to change the energy 
levels so as to stay in resonance with the laser. In two sep¬ 
arate papers (Ref 9 and 10) Phillips and Metcalf, along with 
Prodan, reported results of the two different methods of cool¬ 
ing. Unfortunately in both cases the laser apparently quit 
interacting with the atomic beam. In both cases the most prob¬ 
able velocity was 1100 m/sec. In one case the atoms could only 
be slowed to ~~ 4 O m/yu, . In the second case the atoms were 
only slowed to ~ 6 0 0 . Consultation between the scien¬ 

tists at National Bureau of Standards and Dr. Richard Cook of 
the Air Force Institute of Technology resulted in a theoretical 

investigation of this problem using a one dimensional Fokker- 

* 

Planck equation. This thesis is the result of that investigation. 
Problem Statement 

Given a longitudinal velocity distribution for an atomic 
beam, and the theory of resonance radiation interacting with 
that beam, predict the time history of the velocity distribution 






























using a numerical solution to the Fokker-Planck equation. A 
careful investigation of the one dimensional model developed 
by Cook (Ref ll) is to be accomplished. The model will be 
checked against known experimental results. If the model shows 
the loss of interaction as found experimentally, then a study 
of ways to avoid the interaction loss will be undertaken. If 
the model does not contain the information, other more complex 
models will be briefly reviewed to find avenues of further 
research. 

Current Knowledge 

The background section gives a historical development of 
the important papers in this field. Hence, this section will 
be constrained to dealing with the body of knowledge needed 
for the thesis. This consists of three basic areas. First, 
the initial velocity distribution for the atoms are needed. 
This information is developed in Ramsey’s book (Ref 12). 
Second, a general theory of the interaction of light with an 
atom is required. A general theory of resonance radiation 
pressure of light is developed from quantum electrodynamics 
by Cook (Ref 13). Finally, the transformation of the partial 
differential equation developed from the theory to a differ¬ 
ence equation is taken from Gerald (Ref 14). 

Approach 

This thesis is broken into four main sections. The first 


section is the introduction. The following three sections are 
the main body of the thesis. A bibliography and an appendix 




























follow the main portion of the work. The following is a brief 
description of the main sections in the thesis. 

Section II develops a set of coupled differential equat¬ 
ions that describe the interaction of the laser with the atomic 
beam. These equations are developed using Ehrenfests theorem 
to find the time rate of change of the expectation value of 
the momentum of the atom. It is shown that these coupled diff¬ 
erential equations reduce to a single Fokker-Planck equation 
in certain circumstances. The Fokker-Planck equation describes 
the velocity distribution function of the atomic beam as it 
evolves in time. The Fokker-Planck equation is the major start¬ 
ing point for this thesis. 

Section III describes the numerical solution of the 
Fokker-Planck equation. A difference equation is developed to 
solve the problem numerically. Stability and convergence cri¬ 
teria of this difference equation are studied. A complete anal¬ 
ysis of a particular cooling method is also studied. The re¬ 
sults derived from the numerical analysis are compared to ex¬ 
perimentally derived results for this particular case. 

Section IV draws conclusions and makes recommendations 
for further work based upon this thesis., A review of the major 
theory and results are presented in brief form. An estimate of 
the validity of these results is given. Recommendations are 
made for further experimental and theoretical work in this area 

Assumptions and Scope 


Several assumptions are implicit throughout this thesis. 











The assumptions are detailed below, along with their justi¬ 
fications, so the reader can gain a clear understanding of 
the limits of this work. In all cases a minimum number of 
assumptions are made. Hopefully this will keep the theoretical 
results broad enough so they can be applied to other specific 
cases. The theory developed here and elsewhere in the literature 
assumes that the atom under consideration has only two energy 
levels. At the present time there are no general solutions to 
the problem for multi-level atoms. Experimentally, a magnetic 
field can be applied to the atomic beam to allow only two energy 
levels to be occupied by the electrons (Ref 15). Thus, the 
theoretical results can be compared to experimental results 
without correcting for this assumption. 

Only a one dimensional model of the interaction is con¬ 
sidered. There are indications that this assumption will cause 
some results to be ignored in the model that are present in 
experimental work (Ref 9 pg 1151)• This assumption is made 
since this is the first time that numerical solutions of the 
Pokker-Planck equation for this situation have been attempted. 

An inclusion of two and three dimensional effects cam be added 
at a later time to the difference equation. The major result 
of this assumption is, the finite extent and Gaussian nature 
of the laser beam is ignored. Since the number density in the 
atomic beam is very low, approximately \0** crK* , particle- 
particle interactions in the beam can safely be ignored for a 
first order answer. 

































To allow the Fokker-Planck equation to bo used as 
the master equation governing the atomic motion, it is 
assumed that the external field changes amplitude slowly 
in time compared to 1- ^ . Here f\ is the Einstein spon¬ 
taneous emission coefficient. Since the case considered is 
a traveling wave monochromatic field interacting with the 
atomic beam, the above assumption is justified. If the 
external field is a standing wave, and the atoms are moving 
with respect to the nodal points, the above assumption would 
be unwarrented. 

The scope of this thesis is intentionally limited due 
to time constraints. First the theory of resonance radiation 
is reviewed in the literature and applied to a specific 
cooling case. Second this specific case is analyzed in 
detail. Third the theoretical results are compared to exper¬ 
imental results, and are reviewed for possible future options 
for continuing work in this important area. 































Introduction 


This section contains physical explanation of the 
interaction of a laser beam in resonance or near resonate con¬ 
ditions with an atomic beam. A qualitative explanation will 
be given first to provide a general idea of the physics in¬ 
volved in cooling a beam. A quantitative explanation of the 
physics will follow. The specific explanation is a compilation 
of several articles by Dr. R.J. Cook (Ref 11 and 13). 


General Explanation 

When light interacts with an atom it excites the elec¬ 
trons in the outer shell of the atom to a higher state of 
energy. In 1917 Einstein (Ref 1) gave the first explanation 
of this process. Photons have definite momentum and energy 
described by the famous relations. 


P - Is k Es +\ \w 


where k has magnitude and points in the direction of prop¬ 
agation of the light wave. A photon striking an atom excites 
the electrons when the energy carried by the photon equals 
the energy gap between the electronic transitions. When this 
happens the photon imparts momentum to the atom. After a time 
where “p^ , ft is the Einstein spontaneous emission coeffi¬ 

cient, the atom will decay back to the ground state. At this 
point it can be excited by another photon and receive another 
"kick" in momentum. If the photons are collimated as in a 









laser, the imparted momentum will all be in the same direction. 
Einstein showed that spontaneous emission has equal probabil¬ 
ity of emitting in the positive or negative direction for a 
one dimensional case (Ref 1 pg 67 ). The emission pattern of 
the atom while not isotopic, is symmetric. Thus, the spontan¬ 
eous emission of the atom imparts no net momentum change to 
the atom. A large number of photon-atom collisions will then 
cause the beam to be slowed down if the atoms are propagating 
in the opposite direction of the laser. 

Specific Explanation 

Consider an atom with two allowable energy states, E 
and E a . This atom is exposed to a traveling wave tuned to 
the transition from E, to E a . This traveling wave has an 
associated energy density W. Let P, and P* be the probability 
that energy states E, or E^ is occupied. Einstein proposed 
that in this case the rate of stimulated emission is 

=BWP a . The rate of spontaneous emission is 
R =AWP 4 . Here B and A are the Einstein B and A 

coefficients for stimulated and spontaneous emission respec¬ 
tively. For an atom in steady state condition, i.e. the field 
has been interacting with the atom on 'a time scale long com¬ 
pared to T- a , the rate of exciting processes equal the 
rate of deexciting process. 

Hence, 








The atom receives -*k>k from the absorption and -"kk from 
the stimulated emission. So the time rate of change of the 
momentum is given by equation 1. 

f= -tt * if * - r (d 

This derivation of the force is accurate only for cases 
where the momenta transferred has definite value isk . 

If the field has a variety of component wavelengths, has k's 
pointing in different directions or is not a traveling wave, 
a more gneral explanation is required to predict the observed 
experimental results. 

Two papers by Cook (Ref 11 and 13) provide a more general 
theory of the equations of motion of the atomic beam. The first 
paper treats the motion of the center of mass of the particle 
using Ehrenfests theorem and the optical Bloch equations. It 
neglects fluctuations in the force of the photons on the atoms 
due to the spontaneous and stimulated emission of the atoms. 

The second paper develops a general theory for a two level 
atom in a resonate or near resonate field of arbitrary ampli¬ 
tude and phase. The following will be development of the apr 
proach used for Ehrenfests theorem and the optical Bloch equa¬ 
tions. 

Following Schiff (Ref 15) consider an atom represented 
by a wave function . To find the "equation of motion" of 
this wave packet take the time derivative of the expectation 
value of the position variable. 
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Substitute for the time derivatives using Schrodingers wave 
equation and perform some algebra to obtain, 

Integrate over the second term by parts twice remembering that 
the wave function vanishes at infinity. 

4<*> • -v'C-uYT) *-c 

Perform the differentiation within the brackets to obtain, 

<1“C - <.!>> (2) 

Using the relation, 

<p = d-t . 


In a very similar way one can obtain, 



where V is the potential energy of the system. This is a good 
representation in the macroscopic limit where the finite size 
and internal structure of the particle can be ignored. This 
means the size of the wave packet is small compared to 












































distance over which V changes significantly. Equation 1 and 
2 provide an analogy to the classical equations of motion. 

They also agree with the Heisenberg representation which 
states, 

. x an<) , 

This is because -,£•£ is the potential energy for an atom in 
a field. Here H is the Hamiltonian of the system. For an atom 
in the dipole approximation of an external field, 

fi. 4 A. - 

Here is the electric dipole moment operator and E is the 

— a 

electric field evaluated at a point ^ . H 0 is the Hamiltonian 
for the internal motion of the atom. 

Ehrenfests theorem is equation 4. It is obtained by com¬ 
bining the expectation values of the Heisenberg equations and 
setting, 

F, « ilil> (4) 

Equation 4 is general as long as the dipole approximation 
is accepted. To ease calculations, assume that the external 
field has the form, 































where G is a polarization vector independent of r and *i . 

This form of E is substituted into equation 4 to obtain 
equation 5. 

f = (5) 

This uses the assumption that the wave packet is small compared 
to the changes in the field. This allows the V* £(f.O to 
be pulled out of the expectation value. Therefore in the 
dipole approximation the motion of the atom is due to the 
electric field vector. To find the force acting upon the atom, 
and hence its quantum mechanical ‘equation of motion", consider 
a two level atom with energy states E and E . An arbitrary 
monochromatic field is applied to the system. The form of this 
field is equation 6. 




( 6 ) 































The Optical Bloch Equations 


Several authors (Ref 11 and 16) show that using equation 
6 as the applied field in equation 5, yields a coupled set of 
differential equations. These equations explr.in the motion of 
an atom in an external field. Equations 8 thru 10 are the 
optical Bloch equations. These equations govern the internal 
workings of these atoms. 


^ - rr> V- = ^ “K (_ U V v/L "t V J\j V © ) 

u * (b* e)\>- 4 fly 


v: - i flu 
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Here is the phase of the external field, 

is the amplitude of that field, 


y\- joJL is the on resonance Rabi flopping frequency, 
the detuning frequency, 


ft g ia| t I * is the Einstein spontaneous emission coefficient. 

1*0 ' * ' 

The Rabi frequency is a measure of the field strength at the 
point under consideration. It also represents the frequency at 
which the atom will absorb and undergo stimulated emission in 
a strong field. Equations 7 thru 10 show that the force is not 
an explicit function of position, velocity, or external field. 








































but rather is governed by the coupled set of equations 


shown. The terms involving the Einstein coefficient are 
called relaxation terms. They are obtained thru quantum elec¬ 
trodynamics and are not derived in this thesis. Louisell (Ref 17) 
presents a first principles argument for these relaxation 
terms. 

If the spontaneous terms were not present, the atom 
would oscillate between the upper and lower state at the Rabi 
frequency. The spontaneous emission process, which gives rise 
to the relaxation terms, causes a steady state population dis¬ 
tribution. The physical significance of the variables U,V,W 
are summarized from Allen and Eberly's book (Ref 11). W is 
the single atom population inversion, V is the absorptive com¬ 
ponent of the dipole moment, and U is the dispersive component 
of the dipole moment. 

The Fokker-Planck Equation 

Cook (Ref 13 pg 1087) develops a set of equations similar 
in nature to equations 7 thru 10. The equations developed by 
Cook are more general in nature and are derived from quantum 
electrodynamics. The optical Bloch equations are developed 
from the Ehrenfests theorem and only predict the equation of 
motion of the centroid of the atomic wave packet. The optical 
Bloch equations are consistent with the more general set of 
equations. It is this fact that allows use of the optical Bloch 
equations in the development of the theory. The optical Bloch 
equations represent a less general but still useful starting 



































point for most practical applications. 


An easier way of treating the motion of an atomic beam 
interacting with a laser is to use the Fokker-Planck equation. 
Using the general set of equations developed by Cook, it can 
be shown that the Fokker-Planck equation is accurate to second 
order in (Ref 13 Pg 1096). An explicit derivation of the 
Fokker-Planck equation is given in the Appendix of Cooks paper 
(Ref 13 Pg 1096). The form of the Fokker-Planck equation is 
given in equation 11. 








( 11 ) 


where $ is the velocity distribution function for the atomic beam 
F is the interaction force of the laser with the beam, 
rr> is the mass of the atom and 
D is a diffusion coefficient. 

Since the derivation of the Fokker-Planck equation is 
given explicitly by Cook, it will not be repeated here. The 
restrictions imposed by using the equation, and the development 
of explicit expressions included in the equation are presented 
below. 

The major constraint imposed by the use of the Fokker- 
Planck equation is the smooth field approximation discussed 
by Cook (Ref 13 Pg 1096). The smooth field approximation re¬ 
quires that the amplitude and phase of the external field 
remain nearly constant during a time equal to the relaxation 
time of the atom.This relaxation time is on the order of t-A. , 


















where ft is the Einstein coefficient. If the amplitude and 
phase of the external field are not constant in space, then 
the velocity of the atom must he such that it will not move 
into a region of significantly different amplitude or phase 
in a time "t-’l- . For instance, in a standing wave of wave¬ 

length \ , the atom should not move a distance greater than 
-^during the time X . 

The interaction force in equation 11 is derived from 

equations 7 thru 10. In the case of the smooth field approx 

imation since the phase and amplitude vary slowly in time, 

the values of U,V,W assume steady state values obtained by 
« • • 

setting U,V,W = 0 in equations 8 thru 10. After elimination 
of W one obtains equations 12 thru 14. 

? -- i*»(u<7va t vjx. ve) (; 




(13) 


V = - 3 RJl/ l ©V 4 ft' v 


(14) 


Solving these three equations for F results in, 

^ PI SC ^7 0 •» ^ ( &+ ©) V J\ x 3 

lH(^4©3'+ ft x + 


For a plane running wave of the form, 
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fe=-k-x ^ e.-fc-v ^©--k so j 

f = ft JVWl _ 

Cm (h-H v^ x 4 fV-oA‘3 

This is the force associated with the absorption of momenta 
from photons discussed in the General Explanation section. 
This force is commonly called the spontaneous emission force 
since it is the symmetric emission pattern of spontaneous 
emission, that allows the momenta transferred from absorption 
to effect the atom. 

For a standing wave of the form, 


£lx.O - CCO JY- ^ 6M ftefe ^ 18 

* 

so that, 

f - - ^ ftV SL _ (19 

% 

This force is commonly called the dipole interaction 
force. This force arises from the dipole moment induced by 
the external field in the atom interacting with the gradient 
of the amplitude of the external field. In general, all one 
needs is the form of the external field. This combined with 
equation 15 gives the force of interaction between the light 

































and the atomic beam. 

The diffusion coefficient in equation 11 arises because 
of fluctuations in the interaction force. The diffusion coef¬ 
ficient has two components, one due to the spontaneous emission 
process, and one due to the induced processes of absorption 
and stimulated emission. The derivation of the form of the 
diffusion coefficients is given by Cook (Ref 13 pg 1097). The 
results of this derivation are equations 20 and 21 . 


s H fL£k k—sLi- f 


( 20 ) 


P fxl "fc 4 *** J(wUuctO 

where 

vU''' 5 

Ho b-* 1-.S& 

i- 

C» - a A' 

As shown, both the diffusion coefficients are tensors. The 

. • 

matrix 4 is a diagonal matrix when the dipole moment is 
directed perpendicular to the direction of propagation of the 
external field. For this case and a one dimensional analysis 
\ (Ref 13 pg 1087). If the form of the external field 
is £(«>V* *• «••>( then equations 20 and 21 become equa¬ 

tions 22 and 23 . 
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The physical explanation for the spontaneous diffusion 
coefficient is fairly straightforward. The atom spontaneously 
emits a photon after "being excited in a time period 1* ^ . 
This emitted photon imparts momentum to the atom. Since the 
direction the photon is emitted in is random, many such emis¬ 
sions will give rise to a walk phenomenon. Thus, a random walk 
in momentum space occurs. Following Rief (Ref 18 pg 488), one 
can calculate a diffusion coefficient which is within an order 
of magnitude of that obtained from equation 22. The step sizr 
is V\k » and "the rate ft • Hence, 


Ds- iP'A 




This is within a factor of 5 of that derived from equation 20. 
Equation 20 evaluated in a strong field case yields, 


^1*0'P> 

Since the emitted radiation pattern is dipole in nature 
(Ref 13 pg 1089), the diffusion coefficient due to spontaneous 
emission has a definite tensorial nature. That is, it has 
different values for different directions in space. 

The induced diffusion coefficient can not be so 








































conveniently explained. Absorption and stimulated emission 
processes are not statistically independent. To absorb a photon 
the atom must be in a ground state. To emit a photon the atom 
must be in an excited state. So the probability of whether a 
photon emitts or absorbs is dependent on what has happened in 
the history of the atom. Thus, the statistics of photon absorp¬ 
tion and emission are not Poisonian (Ref 13 Pg 1090). 

The idea that the induced diffusion coefficient can be 
derived from a random walk argument is incorrect in the general 
case. A step size of and a rate J\. leads to an answer that 

disagrees with that found from equation 21 by a factor of £ . 
This factor can in principle be as large as desired since it 
is linear with the field strength. In general the induced 
diffusion coefficient has not yet been physically explained. 

It has been pointed out that a physical explanation can 
be provided for certain very restricted cases (Ref 13 pg 1089). 
For the case of a standing wave external field that is exactly 
on resonance b-o , a physical explanation of the induced 
diffusion coefficient is developed from the optical Stern- 
Gerlach effect. 

When the external field is a standing wave and the laser 
is on resonance with the atomic beam, then there exists two 
separate distributions for the atomic beam. These two distri¬ 
butions are a superposition of pure quantum states for the atoms. 
These two distributions experience a force F 4 and F_ . The 
states are labled , and 5. for the forces they experience. 

The nature of the force is most easily explained by traditional 

















electromagnetic theory. When an electric dipole is in an ex¬ 
ternal field it experiences a force . This is 

the magnitude of the force the atom experiences. The two 
distributions experience forces in the opposite direction. 

From the view of stimulated emission and absorption, the force 
arises from many coherent emissions. The probability that the 
atom is in either distribution is £ (Ref 13 pg 1084). The 

U 

force continues to interact with the beam for a time 
Since the probability the atom is in its upper state is £ for 
the strong field case, the rate at which the atom emitts is 
A/4. This is because the rate an atom emitts is R = AP, where 
P is the probability the atom can emit (£•£=*). 

Thus, when spontaneous emission returns the atom to its 
original level, the entire process can start over. The atom 
can once again feel a force IF. One can see that the beam will 
be split in half each time * . This represents a diffusion 
in momentum space that can be represented as a random walk 
process. 

For a step size in momentum space = and a rate of 

xj , one obtains equation 24. 

3.0"(fcpF r ^ • (24) 

Thus, the diffusion coefficient as calculated from the random 
walk approach is 


0 - * 






( 25 ) 




























Equation 21 gives the induced diffusion coefficient in 
general. For a general standing wave of the form equation 26, 
the induced coefficient is equation 27 . For the case of exact 
resonance c. , one then obtains a coefficient that is de¬ 
rived from the same conditions as before. This coefficient is 
equation 28. 




(26) 


(Wl <27> 


t) c« $ or £= o (28: 

* R 

This coefficient agrees exactly with that of equation 25 
which was developed from the random walk process due to the 
optical Stern-Gerlach effect. It should be remembered that 
for off resonance cases, physical explanations do not yet 
yield a correct induced diffusion coefficient. Equation 21 
is the general form of the induced diffusion coefficient. 

In conclusion, the spontaneous diffusion coefficient 
(equation 20 ), can be explained as a random walk process in 
momentum space size -M* and rate A. The induced diffusion 
coefficient (equation 21 ), does not have a physical explan¬ 
ation that holds for all cases. For on resonance the Stern- 
Gerlach effect can be used to develop a diffusion coefficient 
from the random walk point of view. The total diffusion 


























coefficient appearing in the Fokker-Planck equation (equation 11), 
is the sum of the spontaneous and induced coefficients. 


bi. 


(29) 


Discussion 

As noted above, the Fokker-Planck equation is a special 
ized case for a general set of coupled partial differential 
equations. These equations are presented below in equations 
30 thru 33 for completeness. 
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where Q v ' « a ^ ^ k 1 




and *Uk) x w <5*(u> v ( I - 

'in the dipole approximation. 

Here & is the angle between \ and , 

'rvx is the mass of the particle and 
is the dipole moment. 

Equations 30 thru 33 represent the generalized quasiclassical 










limit for the equation of motion for an atom in resonate con¬ 
ditions. 

It is obvious that the Fokker-Planck equation presented 
in the preceeding section is easier to solve than equation 30 
thru 33. The decision of which equation is to he used is based 
upon the exact problem under investigation. For example, a 
plane monochromatic running wave obviously fullfills the re¬ 
quirements for the smooth field approximation. Thus, the Fokker- 
Planck equation could be used to investigate this problem as 
long as the time scale under consideration was long compared 
to 1/A. On the other hand, a problem with a standing wave 
external field in which the atom moves a significant distance 
compared to the distance over which the field changes amplitude, 
would require that equations 30 thru 33 be solved. So the 
optical Stern-Gerlach effect can not be predicted using the 
Fokker-Planck equation. 

In conclusion, there exists a set of generalized equa¬ 
tions of motion for a two level atom in an external field 
(Ref 13 Pg 1082). These integral equations include all effects 
of field quantization. The quasiclassical limit (Ref 13 pg IO83) 
of these equations are presented here as equation 30 thru 33. 

The use of Ehrenfests theorem and a form for the applied field, 
yield the optical-Bloch equations (equations 8 thru 10). The 
optical Bloch equations are consistant with the generalized 
set of equations but only predict the equation of motion of 
the centroid of the atomic wave packet (Ref 11 pg 227). That 
is, they will not predict such results as the optical Stern- 



























Gerlach effect which is quantum mechanical in nature. For 
the smooth field approximation (Ref 13 pg 1096), equations 30 
thru 33 reduce to the Fokker-Planck equation (equation 11) 
with force and diffusion coefficients given by equations 15 
and 29 respectively. The Fokker-Planck equation is accurate 
to second order in . For the case of a plane monochromatic 
traveling wave interacting with an atomic beam, the Fokker- 
Planck equation is sufficient to solve the equation of motion 
of the atomic beam distribution function. The analysis section 
takes this example and presents the numerical solution of the 
Fokker-Planck equation for this case. 



































III. Application of the Theory to Cooling an Atomic Beam 


Introduction 

This section deals with the application of the theory 
developed in section II, to the problem of using a laser in 
resonance with the atomic beam to decrease its speed. Following 
this introduction, the Fokker-Planck equation will be converted 
from a partial differential equation into a difference equation. 
The numerical solution of this difference equation is studied 
in view of the constraints imposed upon the step sizes of 
time and velocity, due to convergence and stability of the 
numerical solution. A description of the particular case studied 
and its initial conditions are then given. The difference equa¬ 
tion is solved numerically and compared to experimental results. 

The Fokker-Planck Equation as a Difference Equation 

As shown in section II, the Fokker-Planck equation serves 
as a method of calculating the time history of the atomic beam 
distribution function in cases where the smooth field approx¬ 
imation applies. One method of numerically solving this equation 
is to convert it to a difference equation and solve using the 
method of finite differences. 

If the differentials in the Fokker-Planck equation are 
replaced by difference quotients, the resulting equation is 
known as a difference equation. Deriving these difference 
quotients is straightforward. If there exists a function “S 
that has continuous fourth derivatives, then it is expandable 
in a Taylor series expansion as in equations 3^ and 35. 
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Here K is the step size. To find difference quotients for Vtx*') 
or , one solves equations 3k and 35 for the derivative 

wanted. For instance, if equation 35 is added to 34 and 
subtracted from the resulting sum, one obtains equation 36 . 

* UvrO - s S‘*u*w ^ 4 (36) 

Here %(x n *W'> means the function evaluated at the value , 

where K is an arbitrary step size. This is rewritten in equa¬ 
tion 3? as, 

- ^H*^4 S(x,>-V) ^ (37) 

V ^ 

O(V) means that the error approaches proportionally to W* 
as V>~* O . This is known as a central difference formula for 
the second derivative (Ref 14 pg 219)» since the derivative 
at point 'X* depends on values of the function at velocity 
points before and after it. A forward difference formula for 
the first derivative is derived from equation 3k by solving 
directly for as in equation 38 * 
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As shown, this has errors proportional to K . Finally, a 
central difference formula for the first derivative is given 
in 39. 


S £(x,y»V>V + (39) 

2 K 

Corresponding to equations 37 thru 39» are similar equations 
with the variable "t replacing x . Obviously, the choice of 
variable is immaterial. 

To allow use in a partial differential equation, the 
difference equations must be modified. To solve the Fokker- 
Planck equation in one dimension, assume arbitrary step sizes 
for time and position. Thus, $(***0 becomes fcXn') , 

and S(itVT) becomes . Partial derivatives with 

respect to time and space can then be written as equation 4-0 
thru 42. 
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, forward difference. (40) 

, central difference. (41) 

, central difference. (42) 


Again the choice of variable is of no substantial consequence. 
The only assumption is that the function has continuous fourth 
derivatives for the variable in question. 



























Many difference equations can be written for one partial 
differential equation by using different combinations of for¬ 
ward, central, and even backwards difference equations for 
the partials. The Fokker-Planck equation must be solved for 
successive increments in time and velocity. Thus, a forward 
or central difference must be used. To ease calculations and 
file space, a forward difference will be used for the partial 
with respect to time. The second derivative must be written 
as a central difference if one wishes to avoid using a first 
derivative to calculate the second derivative. The choice of 
forward or central difference for the first partial with re¬ 
spect to velocity is at this time arbitrary and will be left 
for later. 

A time and velocity plane is constructed to solve the 
Fokker-Planck equation. The axes are scaled from 0 to as large 
as desired in increments of and . The partial de¬ 

rivatives presented in equations 40 thru 42 are written again 
in equations 43 to 45. 
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The superscripts denote the time position at which the function 





















The Fokker-Planck equation is a parabolic differential 
equation (Ref 19)• Thus, the equation is solvable using the 
explicit method of finite differences (Ref 14 pg 220). This 
means that when the partials are replaced by the difference 
quotients, the value of the function at the ith velocity can 
be solved explicitly for the j + 1th time. This is evident 
from 46 and 47. All values in the equation are known except 
for $?' . 

The initial atomic distribution function is given by 
Ramsey (Ref 12). This represents the atomic velocity distri¬ 
bution funtion as it comes out of the atomic oven. This 


















distribution function will be changed as time proceeds as 
prescribed by equations 46 and 4?. The general forms for 
and D? are given by equations 15 and 29. The explicit values 
used for the force and diffusion coefficient depend upon the 
exact problem to be solved. The next subsection discusses the 
particular case to be examined and derives explicit equations 
for F* and fo; 

Slowing an Atomic Beam Using A Scanning Frequency Laser 


Consider a monochromatic plane traveling wave counter- 
propagated against a collimated atomic beam. The form of the 
wave is given in 48. 
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The force of interaction between the field and the beam is 
given in 49. 
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The total diffusion coefficient was derived in equations 22 
and 23 . It is given here as equation 50. 





















The wave vector k«- IklU , since the laser is counter- 
propagated against the beam. 

In both equations there is a term kv . This is the 
effective detuning since and - kv is the doppler 

shift term. The laser starts at a frequency tuned to resonance 
with the initial velocity atoms V* . The effective detuning 
is zero at the peak of the force. As the laser frequency is 
swept, different velocities come into resonance with the laser. 
That is, different velocities make the effective detuning equal 
to zero. The force is maximized when the effective detuning is 
zero. 

k-vk\rpzO , ATp represents the velocity at peak force. ( 51 ) 
At time t* O , the velocity v:v„ > so that 

- Vat 0 (52) 


The force then assumes the form 


T - V 

+ 3JV 


(53) 


If the laser frequency changes, then & changes in a 
manner yet to be determined. As & changes, different velocities 
will satisfy equation 51- Thus, the force will be resonant 
(i.e. a maximum), for different velocities of the atomic velocity 










distribution function as the frequency of the laser is changed. 
The force is opposite the beam so the atoms are slowed. If the 
laser frequency is changed so that the atoms in the beam slow 
just enough to remain at the peak of the force, then 

p- m> \r - - U k (54) 

Integrating 5^ yields a velocity at which the force is maximized 
at 


"irp - ft Jl* kt (55) 

Here v* is a constant of integration and equals the initial 
velocity that is on resonance. Equation 51 gives another relation 
for the velocity at the peak of the force. Combining these two 
yields 

- kv. - - V W, - (56) 

So if the atoms stay in resonance with the force, then the 

detuning changes with the scanning frequency of the laser as 

» 

k = - kv t -v Kt , ^ , R A x k v (57) 


This is the detuning rate at which the atoms slow just enough 
to stay in resonance with the laser. If the rate is any faster 






































the atoms will not be slowed enough to stay with the peak 
force. The detuning rate should be less than the maximum cal 
culated in 57. This indicates the inequality 


R ^^ . R has units of s«.c in MKS. (58) 

The physics of the problem are now established. A laser 
interacting with an atomic beam which is characterized by its 
distribution function % , exerts a force given by 49. The 
laser initially is in resonance with an initial velocity tr* . 
The laser frequency is swept at a rate given by 58. This causes 
atoms to slow down and the velocity distribution function is 
changed. A diffusion coefficient given by equation 50. causes 
heating and moves atoms out of the peak of the force. As the 
laser frequency decreases, substantial numbers of atoms are 
slowed to lower velocities. The final velocity the atoms are 
slowed to is proportional to the doppler shift kv* 

The force in equation 49 is Lorentzian in w . Thus, it has 
a characteristic full width at half maximum spread in frequency 
of iR'OjlV*- . This corresponds to natural and power broad¬ 
ening of the force due to the atomic response (Ref 11 pg 226). 
For typical numbers, this is much narrower than the FWKM of 
the initial distribution function. So as the laser is swept 
in frequency, the distribution function becomes narrower and 
moves to a slower velocity. 





























Stability and Convergence Analysis 


As mentioned before, the step sizes for ax and 
are arbitrary within limits set by stability and convergence 
criteria. The convergence criteria are derived using the prin¬ 
cipal that the numerical solution approaches the real solution 
as k'v and tend toward 0. Stability is defined as the 

damping of error terms as they propagate through the calcul¬ 
ations in time. The two difference equations introduced in 
equations 46 and 47 are investigated for convergence criteria 
first, then for stability. The two equations are then compared 
and the one that has the least stringent requirements on 
and is chosen to solve the Fokker-Planck equation. 

The equations needed for the convergence and stability 
analysis are equations 59 thru 61. Equation 59 is the Fokker- 
Planck equation. Equation 6 0 is the Fokker-Planck equation 
rewritten as a difference equation, with a forward difference 
used on the first partials, and a central difference used on 
the second partial with respect to velocity, henceforth called 
method 1. Equation 6l is a difference equation, with the ve¬ 
locity partials written as central differences, and the partial 
with respect to time, written as a forward difference, hence¬ 
forth called method 2. 
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As before, the superscripts refer to time increments and the 
subscripts refer to velocity increments. From now on, the 
convention will be used that all superscripts will be dropped 
unless they do not equal j. Therefore, all difference equations 
will have only one term with superscripts. The convergence 
analysis for methods 1 and 2 follow. 

The exact solution to the differential equation is *Un 
and the numerical solution is 4 a • The error involved with 
substituting the numerical for the real solution is 


Thus, 




, and 


To analyze the convergence constraint for method 1, substitute 
the above relation into 60 and collect terms to obtain 
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Here <•* .A Taylor series expansion of •uL, *Ua‘ 
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Substitute the above relations into 62 and cancel the appro¬ 
priate terms to obtain 63 . 
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Here ^( v « ^ ^ [- &. * |^*v (£„0*£ 
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with all r* expressed as r* 

The magnitude of the maximum error in time t 1 * 

is E 1 *'- T J . To calculate this error, take the maximum 
values in equation 63 and the magnitude of ^ . 

‘"•[Mfr E “ # E ^ ^ riv]. E‘[ 1-igJJ* ,£^1 (6 
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where t)!*» and are the maximum values of ft* and F« 

This leads directly to 65 . 
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If the coefficients of E 4 and |ftl are positive, then 
the following inequality holds. 

E 4 silnl (66) 

and thus, k.'tlP'N) . Applying this result at 

different times leads to 6?. 

E^'< £^*+ 3 ut\iA\ . i bilhl (67) 

The last term kMlM- since E^iCi 

at time E * = ^ . 

Equation 67 shows the difference equation will converge 
to the real solution after j iterations. This will happen 
when initial errors are zero and when is arbitrarily small 

with lI'M finite. The initial condition in this problem is the 
velocity distribution function, which is theoretically and 
experimentally well known. So the initial errors can be made 
zero. The term can go to 0 as goes to 0 as long as ItM 

is finite. is given in 63 and by inspection it is finite 

as long as and *v% 4 . 

As pointed out above, the coefficients of E 4 and lI'M 
must be positive for 66 and 67 to hold true. So the convergence 
criteria are found from requiring that the coefficients be 
positive. Therefore, 



























Now Dmax ■ = jfe** [*'»■* [>*' *1 IS) 


and Fmax = F m t "L. 
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Obviously Dmax can be made positive since the Rabi frequency 
can be made larger than the Einstein A coefficient. Equally 
obvious, the term Fmax is always negative for this case. Thus, 
the first and third equations can provide no constraints upon 
the difference equation. The second equation is the only equation 
that can provide such a constraint. To find this constraint, 
solve the inequality for r. Replace r with its value . 

Upon doing this one obtains 




These are the convergence criteria for method 1. 

% 

The convergence criteria for method 2 are found in much 
the same way. Let the real solution corresponding to the differ 
ential equation be . The solution to the difference equa¬ 
tion is . The error involved is . This leads again 

to h * *\>1- 0* . 






















Substitute this into 6l and obtain 69. 


e* = Ain +1 U ■ V 


Expand the terms Ml, » , as Taylor series as 

before. Collect all terms in and express them as M. Equation 
69 then becomes 70 . 
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+ i .r.gj-,Y| + ht t“\ (70) 
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The maximum magnitude of the error for the calculation at time 
V is E' , so that once again, 

e‘* ^iia\i E i_, o*Uh\i • •• < *t\v\\ 


This is true if all coefficients of the error term are positive. 
This gives the constraint upon r such, that the solution from 
the difference solution converge to the solutions of the differ¬ 
ential equation. 
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Equations 68 and 71 represent the constraints upon methods 1 
and 2 respectively for the solutions of the difference equations 
to converge to the solutions of the differential equation. 

The stability constraints are now obtained. Expand the 
velocity cistribution function as a fourier lattice series 
(Ref 20 pg 3)• One then obtains 72. 

^ » where "U„(UV ^ l ®^ n (?2) 
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The subscripts denoting velocity have been changed to n since 
a complex term has been introduced. Equation 72 shows that 
the distribution function at velocity n is a series of fourier 
components. Each component in this expansion is a solution to 
60 and 61. If the coefficients are constant then each component 
of this series can be seen to be a eigenfuntion of equations 
60 and 6l. Substitute the component into 60 and obtain 

l uwf. . UsdaS"] (73) 
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Obviously, the term in brackets is the eigenvalue. So each 
component of the expansion conforms to equation 74. 

• 'X'Oo (?4) 

Equation 60 and 61 show that each value of the velocity 
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distribution function for time tit** is calculated using 
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existing values of the distribution function at different 
velocities. If an error is introduced as calculations proceed, 
because of round off errors or other problems, then it affects 
subsequent calculations. The error is 

, where is the value that is incor- (75) 

rect at time j and velocity n. This error term will give 
solutions to equation 60 of the same form as 7U except that 
*\An V ' will be . Starting with 7^ and instead of 

'M?' at time 

For the next time iteration one obtains 

( 76 ) 

Here the superscripts on the e*s are time increments and the 
superscripts on the X's are exponents. For each time iteration, 
the error is multiplied by the eigenvalue X . 

Stability for numerical calculations means that as succes¬ 
sive calculations are made, the errors tend to damp out. Looking 
at 76, one sees that after j iterations 

* (77) 

If due to the initial error , one requires that after j 
iterations Ot Vt 4 , then obviously lX\ < \ . If W\ > ' 
then the errors grow instead of decrease. To find the requirement 
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for stability then, one requires that 


r •»". 


The constraint upon method 1 is then 


| 1 * 1 '®| 4 1 (78) 

This simply states that the eigenvalue found in 73 has a 
magnitude less than 1. Euler's relations have been used when 
writing 78 to turn the exponents into sin's and cos's. 

Calculating the indicated magnitude is tedious but 
straightforward. The real and imaginary parts of the equation 
are squared and added together. The square root of the resulting 
quadratic is the magnitude. After this, one solves for r and 
obtains 79. 
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This calculation is made for one particular . 

Since the value of the velocity distribution is a sum of a 
series of the s , each must remain stable. The errors 

must not be growing in any of the components of the series, 
or these errors will swamp out all the other components. For 
every value of § , equation 79 must be true or the difference 
equation will be unstable. Thus, every value of © that gives 
a unique value to 79 must be satisfied. Following Kittel 
(Ref 20 pg 3)» 'the allowed values of 0 are given by the 
periodic boundary condition that 





















So 
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This will only he true if or t\ . 

So the values of © that are allowed are 

-t\ < © <■ ♦ n or -T»<kvVt^ - n ( 81 ) 



Since 79 must he satisfied for all values of ^ , only 
the most stringent requirement on the inequality is important. 
Since «■* in 79. it can he seen that what this equation 

really means is that *t must he smaller them a definite 
number for every value of fv . The most stringent requirement 
corresponts to the smallest zA that is calculable. To min¬ 
imize the . the (1-cos © ) term must he maximized. The 

maximum value for this term is 2. Place this in 79 to obtain 
the constraint on ft for method 1. 
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which becomes ft < i e>vmTfQ} upon cancellation of 

appropriate terms. The terms cancel since Fmax =- Wmaxl . 






























Inspection of 68 and 82 show the criteria for convergence 
and stability are the same for method 1. 

The stability criteria for method 2 are found in a 
similar manner. Insert equation 72 into 6l to obtain 
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( 83 ) 


Using Euler's relation 

c\* 1 ^ - yr ^\j-FrO»A^ 1 , where ^vV» (84) 
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The term in brackets is the eigenvalue. As before its magnitude 
must be less than 1 for stability. So the constraint upon 
method 2 is 

(i V -|\ < 1 1. w e>\. < 1 (85) 
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Solving this inequality for r and remembering that r = kU v , 
one obtains 86. This is the constraint upon the time increment 
vs the velocity increment for method 2 due to stability require 
ments. 

H fcv*’ 
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Once again, all terms of the expansion must obey the require¬ 
ment that they not have errors growing as calculations proceed. 
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Thus, all values of cos 0 must be examined and the most 
stringent satisfied in 86. The maximum value for cos © is 
+ 1, and now 


tX 4 . 



for cos © =1, or (8?) 


tX 4 ) 

2 


for cos 0 =-l. 


These criteria agree with Roach (Ref 21 pg 45). 

The above stability analysis is valid only for local 
conditions. An eigenvalue must be constant for a given problem. 
Both the eigenvalues 73 and 84 contain the terms Fmax and Dmax. 
In the original equation these terms depend on velocity and 
time. The constant values of-the peak force and diffusion 
coefficient were used so that the calculations could be per¬ 
formed. Thus, the stability analysis is really only valid 
for regions where the force and diffusion coefficients remain 
constant with respect to their current values. This require¬ 
ment is met if the step sizes for velocity and time are small 
enough to limit the change in force and diffusion. So overall 
constraints on the step sizes are given by the requirement 
that the force and diffusion coefficient can not change dras¬ 
tically over the step size. This leads to the term local 
stability. If the local stability conditions are not met, 
then certainly overall stability can not be insured. Thus, 
local stability is a necessary, but not sufficient condition 





















for overall stability. 

The two methods are now compared to determine which 
has the least stringent requirement on the time vs velocity 
step size. An initial step size will be used, that insures 
the force and diffusion coefficients do not change significantly 
over the step size. The most stringent constraint on each 
method will be compared to see which offers the largest time 
vs velocity step size. Comparing the most stringent constraints 
will insure that all stability and convergence criteria are 
met. Picking the method that has the least stringent of these 
constraints promises faster run times. 

Equation 88 shows the requirement for stability and con¬ 
vergence for method 1. Equations 89 and 90 give the require¬ 
ment for convergence and stability respectively for method 2. 
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Prom before we know that Dm is 
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Fm = 
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For any realistic case, Jl is 5-10 times larger than A. 
After discarding insignificant terms, an order of magnitude 
calculation for on resonance will yield 


= (.■K kV ft and 
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Fm = - Avft x kk c - 
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For comparison purposes of the two methods, values are needed 
for Dm and Fm. Using the following values, one obtains 


Dm = 3-T % > 10' 
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Fm = - S.il % lCf a ' < * n't 

where U = l * lo * 7 rv^ ’ 

k, = ».osM 3-s.« t 

A = l 
m = H 

Methods 1 and 2 can now be directly compared to see which 
has the least stringent requirement on the time increment vs 
the velocity increment. For method 2, the constant value in 


90 is 
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Method 2 ( 95 ) 
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For method 1 the calculation is more tedious, A table of values 
calculated from 88 , of av vs is given in Table I. This 

shows that method 1 has a less stringent requirement upon the 
time increment. So method 1 will be used as the difference 
equation to solve the Fokker-Planck equation for new values 
of the velocity distribution function. 


Analysis of Scanning Frequency Laser Cooling 

The theory of using a scanning frequency laser to slow 
the speed of an atomic beam was developed earlier. A specific 
case is now analyzed using this theory. Prodan and Phillips 
(Ref 10 pg 137) report the use of a scanning frequency laser 
tuned to 589 nm (Ref 22 pg 3826 ) to slow an atomic beam. The 
laser is tuned so as to initially interact with atoms at a 
velocity of 1100 m/sec. The laser is swept at a rate of 6.4E11 
Hz per sec. The laser intensity is 5100 watts/m2 (Ref 23 ). 
Experimental results were reported for scans of 480 and 750 
MHz. This information will be used to generate a case similar 
to this so comparisons can be made between theoretical and 
experimental values. 

Several important constants are needed. The Rabi 
flopping frequency is given in 96 . 
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This is because the dipole moment 
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(Ref 13 Pg 1086), (97) 
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Here, I is the intensity of the laser in watts/m a . Usdng 
Prodan and Phillips numbers then, JV. = S.Snicr sec" . This 
gives a maximum allowable scan rate of iauio' 1 sec"' from 
equation 58. Prodan and Phillips reported a rate of 4.02E12 
sec "'’and 6.28E12 sec**'. Both of these are slower than the 
calculated maximum value. Other needed numbers are 


K * -3^- - \ > \Q 7 rvC ’ 
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( 100 ) 


Table 1 gives acceptable values of the time increment for 
various values of the velocity increment. If the force and 
diffusion coefficients are not to change by 1/10 over the step 
size, then the velocity step must be less than 5 m/sec. A 
step size of 1.5 m/sec for velocity and 52-7 sec for time will 
be used. 

The program used to solve the Fokker-Planck equation 
is given in Appendix I. A brief description of the program 
will be given here. Figure 1 is a flow chart for the program. 
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The program is written in Fortran 77 and was run on the AFIT 
Vax computer. 

The program has three start modes. One exits the program 
immediately. The second starts the program for the first time, 
and the third restarts the program if it was broken off in 
the middle of a run. The program is sometimes terminated 
before completion because the running time becomes excessive 
for various combinations of ai , c-u , and the rate of de¬ 
tuning. When the program is restarted, the data vectors are 
read back into working memory. These vectors are saved as a 
routine process within the program. The program then skips 
directly to the section that solves the Fokker-Planck equation 
using the data vectors and method 1. 

The initial start mode reads in 21 values which charac¬ 
terize the problem. The data vectors are then created. The 
first data vector is the velocity. It is calculated using a 
simple summing process and is called VSLCTR. The initial 
velocity distribution function is calculated using equation 
101 (Ref 12 pg 20). It is named DFTN, 




( 101 ) 


where = 


■ KP -]' 1 


and v= Velocity (m/sec) 

k = Boltzmans constant, 
T = Temperature * K 
m = Mass of atom. 
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The force is computed using equation 49 and is called FORCE. 

The diffusion coefficient is calculated using equation 50 
and is called DIFUSN. 

The section that solves the Fokker-Planck equation using 
difference equation method 1 is the major part of the program. 

The values for the distribution function are calculated for 
the next time increment using values at the current time. 

Thus, the driving mechanism for the program is an iteration 
loop in time with counter T . Instead of recomputing the force 
and diffusion coefficients at each time increment, a interpol¬ 
ation routine is used to shift the values down the velocity 
axis an amount equal to the amount they move in the time 
increment. This can be done because the force and diffusion 
coefficients are represented by a curve dependent on the 
velocity that is translated to lower velocities as time proceeds. 
An examination of equations 49 and 50 will convince the reader 
of this. 

The change in the distribution function called DELDFN, 
is calculated for the ith velocity and the j + 1th time using 
the method 1 difference equation. This change is added to the 
old value of DFTN and a new (i.e. j + 1th) distribution funtion 
is calculated all along the velocity axis. The area under the 
distribution function represents the intergral of the curve 
over all velocity space. Thus, it represents the number of 
atoms under consideration. This number should remain constant 
since the number of particles in the beam is a constant. The 
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value for the area is called TOTNUM and is calculated using 
a rectangle integration formula. 

The results of the calculations are printed into two 
different files. A table of output values called WRTDAT1, is 
used as the input to the graphics package. The table has four 
columns with no headers or other identifying information. The 
four columns are the VSLCTR, DFTN, FORCE and J values respect¬ 
ively. The second output file is called outputx where x is a 
number th?.t keeps track of the output file number. This file 
contains the time increment and the value of TOTNUM, as well 
as the values of velocity, DFTN and FORCE. The outputx and 
WRTDAT1 files contain values for five different time increments. 

It should be noted that the do loop that runs the outputx 
and WRTDAT1 files, prints every Mth velocity value. This cam 
significantly cut down on the file space required. It does 
tend to cause losses in the accuracy of the results since the 
maximum value for the distribution function may be in between 
the Ith and I + Mth values. If K = 1, the stability and con¬ 
vergence criteria suggest that the maximum value will be printed 
for all the values at all the velocities. 

The values for TIME and TOTNUM are printed every 100 
time iterations. A counter called Q is incremented each time 
step. When it reaches 100, the data vectors discussed above 
are stored and the values for TIME and TOTNUM are printed. The 
value of TOTNUM gives a running check of the stability of the 
difference equation. If TOTNUM changes significantly, it is a 








































sure sign that the difference equation is unstable. 

To insure that the convergence and stability require¬ 
ments are met, the program will be rerun with half the step 
sizes for the velocity and time increments. If the values 
of FORCE, DFTN, or velocity where the peak of the distribution 
function sits, change much with respect to their initial 
amounts, then the step size will be rehalved untill it changes 
no more than 20#. The following is a description of the results 
for the difference equation for various initial conditions. 

If the laser is not swept in frequency, then the force 
interacts with only one velocity of atoms. Figure 2 shows the 
initial velocity distribution function for the atoms before 
the laser is turned on. Figure 3 shows the distribution function 
after the laser has been turned on for one millisecond. As 
shown in figure 3» "the force causes a narrow spike to appear 
in the distribution function. The peak is very narrow compared 
to the initial distribution function. The computer took ap¬ 
proximately 950 CPU seconds for this calculation. 

If the peak is Gaussian, then the FWHM is related to 
the v spread of the peak by 102. 

(fvsuiy s a ■ (l02) 

The <t is the root-mean-square deviation of v from the mean -v 
of the Gaussian distribution (Ref 18 pg 24). It therefore 
provides a direct measure of the width of the distribution 
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function about the mean -vr (Ref 18 pg 34). The «■ is important 
because it can be related to a characteristic temperature 
for Gaussian velocity distribution functions. 

Following Reif, the o* is related to this temperature 
by 



Ik "iSU* V 

where k is Boltzmans constant and m is the mass of the parti¬ 
cles (Ref 18 pg 266). So the narrower the peak in the distri¬ 
bution, the "colder” the temperature of the atoms. Using equation 
103, one can compare the initial and final temperature of 
the atomic beam by knowing the FWHM of the distribution function 
before and after the laser interacts with the beam. The initial 
FWHM of the velocity distribution in figure 2 yields a tem¬ 
perature 543.6 °K. The final FWHM of the velocity distribution 
from figure 3 yields a temperature of ,2°K. As can be seen, 
thfe beam is "cooled" by the laser interacting with it. The 
laser soon moves all the atoms in resonance with it to slower 
velocities. If the frequency of the laser does not change, 
the atoms will assume a steady state solution much like that 
shown in figure 3* 

The atoms speed is not slowed appreciably for this case. 
This is because the force function is very narrow. When the 
atoms slow by a small amount, they go out of resonance with 
the laser. To get a feel for how small a change in velocity 














this is» look at the FWHM of the force. The FWHM of the force 
is (Ref 11 pg 226). 

FWHM = (PfO/tV’- = 1.41 *10* sec"' . (104) 

This FWHM has units of radians/sec. In velocity this converts 
to 

*>v = = 14.1 m/sec. (105) 

Thus, the atoms only have to slow ~15 m/sec and they are out 
of resonance with the force. This is why the laser frequency 
must he swept. 

To see the effect of the laser frequency sweeping, a 
case similar to that of Prodan and Phillips will he run. The 
laser is swept for 480 MHz at a rate of .64 GHz/millisec. This 
corresponds to a rate of 4.02E12 sec'*" in radian measure, 
with a total time of 7.5xio“ s sec. Figure 4 and 5 give the 
initial and final positions of the distribution function. 

As can he seen, the laser sweeps the atoms to a slower overall 
velocity as well as narrowing the peak. The peak of the dis¬ 
tribution function is at a velocity of 740 m/sec. The time 
and velocity steps have been cut. The step size was halved 
since the value for the peak of the distribution function 
changed by 27 ft after being halved once. This was felt to be 
large, so the step size was repeatedly halved. This time the 






























































































resulting change in the peak of the distribution function was 
23 fo. The FWHM of the narrow peak is 7.03 m/sec. This corres¬ 
ponds to a temperature of 2.48*iCf x *K. 

Another case experimentally run by Prodan and Phillips 
was with the laser sweeping 750 MHz at a rate of 1 GHz/millisec. 
This corresponds to a rate of 6.28E12 radians/sec x and a time 
of 7sec. Figures 6 and 7 show the results of the dif¬ 
ference equation for this case. Again the overall results are 
the same. The atoms are slowed down and the distribution 
function becomes very peaked. The FWKM is 6.5 m/sec, which 
corresponds to a temperature of 2.122*10"' °K. The final 

peak of the distribution function is at velocity 586.75 m/sec. 

If the atoms were following the force exactly, theory would 
expect the final velocity to be 660 m/sec. 

An interesting calculation is to see if there exists a 
lower limit on the translational velocity of the atoms and 
if there exists a lower temperature bound on that peak. A 
rate of 4.02*\<o' x sec' x was used for a time of 2.475 - ** 0 ’ 1 sec. 
Figures 8 thru 13 show the atoms slowing down as the laser 
frequency sweeps down. The calculations show that the trans¬ 
lational motion of the atoms can be stopped completely. Indeed, 
they can be turned around and be made 'to go the other direction. 
However, a lower limit is found to the temperature. The smallest 
FWHM measured from the calculations was 6.25 m/sec. This 
corresponds to a temperature of 1.96 *K. 
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Prodan and Phillips of the National Bureau of Standards 
have reported experimental results for a similar case to that 
just analyzed. In that experiment, the laser was focused an 
amount equal to the divergence of the atomic beam (Ref 23 ). 

Thus, the electric field is always perpendicular to the atoms 
velocity. This was done so as to make the laser look more like 
a plane wave. Table II sumarizes their results and compares 
them to the results from the analysis performed before. 

Included in Table II are theoretical calculations for the 
velocity at which the peak of the velocity distribution function 
will end. These values are obtained from 106. 

fcvs (106) 

U 

As can be seen from Table II, the values for the end 
velocity of the peak agree with the experimental values 
within 10j5. The values for FWHM do not show such agreement. 

As will be shown later, the guassian nature of the laser 
leads to a different force term from that used in the 
analysis. This different force is analyzed to order of 
magnitude to see if it should cause changes in the final 
value of the FWHM. The conclusion at this point is that 
the guassian nature of the laser does not make significant 
changes unless focussed to a very small spot size. Prodan 
and Phillips do not give experimental error bars for the 
values of the FWHM (Ref 10 pg 140). Thus it is possible 






































































that the factor of six difference between experimental 
and numerical values for the FWKM is not exceedingly 
large. Furthur data is required for a more definite 
resolution to this descrepancy. 

The important consideration is that the numerical 
calculations predict the overall characteristics of the 
experiment quite well. Only one major effect seen exper¬ 
imentally is not shown by the numerical calculations. This 
effect is best seen by looking at figure 14. This is a 
reproduction of Prodan and Phillips data (Ref 10 pg 139). 

As can be seen, the values for the peak of the velocity 
distribution function decrease for a further slowing of the 
beam. A review of figures 8 thru 13 shows this effect is 
not seen by the numerical calculations. This is an impor¬ 
tant consideration since this effect is responsible for 
current lower limits obtainable experimentally for trans¬ 
lational velocity reduction. At this time the descrepancy 
between experiment and numerical calculation is unresolved. 
The numerical calculation predicts that the beam can be 
slowed arbitrarily close to 0 m/sec in translational speed. 

As seen in Table II, the values for the FWHM for 
the two different experimental cases remain nearly constant. 
This is explained by the fact that Prodan and Phillips ran 
the 750 MKz sweep at a slightly faster rate than the 480 MKz 
sweep. This resulted in the interaction time for the laser 
and atoms being equal for the two cases. For both runs, 






















































sec. 


the laser interacted with the atoms for 7.5*l6~ 

This suggests the FWHM is a function of total interaction 
time. To check this hypothesis, another computer calculation 
was performed for the 480 MHz sweep. In this case the rate 
was increased to 6.28•xic/ 1 sec'"*’. As can be seen in figures 
15 and 16, the final velocity of the peak is at the same 
value as before. The FWHM of the peak has changed to 
6.644 m/sec for this case. This is a reduction in the FVJHM 
of 12.1# from the previous case. 

The physical cause of this effect is not yet under¬ 
stood by the author. It does lead to the intresting 
prediction that for a given frequency sweep, a faster 
sweep rate will yield narrower FWHM at the end of the 
interaction. 

Cook gives the lowest achievable temperature for the 
case of a standing wave laser cooling an atomic beam using 
the dipole force to slow the atoms (Ref 24 pg 979). To 
within an order of magnitude estimate, similar arguments 
yield a lowest achievable temperature for the traveling 
wave case of 


j. 


tVv <S X - 


* A 


(107) 


Here the 0“ is related to the temperature by equation 102. 
This yields a smallest expected FWHM of 
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The above formula yields a value of the minimum allowable 
FWHM of about 1.7 m/sec. Table II shows the numerical values 
of FWHM are in agreement with this calculation in that they 
do not violate it. 

As mentioned before the force induced by a guassian 
laser is different from that calculated for an external 
field in the form of a plane wave. To see this difference 
the following is presented. 

Consider the form of the force in equation 15 
evaluated for an external field of a guassian laser. The 
field of the laser is presented in equation 109. 

^ CWt- vwt} (109) 


Here is the beam waist and 

is the field strength, A t >A* ^ 

When this inserted into the force equation, one obtains 
a translational force the same as before. A transverse 
force is also present. For on resonance ( b =0), the 
transverse force has the form of equation 110. 
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If the laser is counter propagated to the atomic beam it 
will defocuss and move the atoms away from the on axis 
position (Ref 11 pg 20). This is because for a counter 
propagated laser . Thus, since the gradient points 

inward, the direction of the force is outward. If this 
force is of the same magnitude as the translational force 
then the on axis velocity distribution function will have 
a smaller value. This could change the value for the FWHM. 
This is a critical question since the numerical values of 
the FWHM are smaller than the experimental values. To 
resolve if the guassian nature of the beam should be 
neglected a comparison of the relative magnitudes of the 
translational and transverse forces is performed. For., 
this: comparison the maximum value of the transverse force 
is needed. This is computed much as before. The only 
difference is that a gradient must be taken of the new 
Rabi frequency. After this is accomplished one obtains, 

The maximum is now calculated noticing that ki > A » r J\ , 

and that a maximum will occur when . This leads 

to a maximum transverse force of, 



Comparison of numerical values with the maximum translational 

- -V 

force shows th$H; U> % = \0 i* for the two forces to be equal. 




































Thus the laser would have to he very focussed for the two 
forces to he of the same magnitude. This shows that for 
lasers with a waist of 1 cm the transverse force is negligible. 
Thus the use of a plane wave for the external field should 
not affect the FWHM value. 

The amount of spread in the peak of the velocity 
distribution function is a measure of its temperature. If 
the velocity and the spread in the velocity is small enough, 
the atoms can he inserted into a trap for neutral atoms. 

This trap is described by Ashkin (Ref 25 ). Cook and Hill 
have a new type of neutral atom trap using the "radiation 
force exerted by the thin evanescent wave that is generated 
on the surface of a dielectric medium when laser light is 
totally reflected internally at that surface"(Ref 26 pg 258 ). 
The requirement on the translational speed for both these 
traps are -*5 m/sec. Using a one dimensional analysis there 
seems to be no reason why these requirements cannot be met. 

A discussion of the amount of confidence warrented 
in the numerical calculations is now in order. The stability 
and convergence analysis govern the step size for time and 
velocity. To meet the smooth field approximation, and 
hence the condition for local stability, the step size of 
velocity was limited to less than 4 m/sec. After the 
calculations are performed, the step sizes are halved to see 
if there are any changes in the answers. When this is 
accomplished, it is found that changes made independently 
































for the time increment, do not effect the results as long 
as the time increment is less than its maximum allowed by- 
equation 96. Changes in the velocity step size are found 
to effect the answers. This indicates that the stability 
of the numerical calculations are assured but that some 
question convergence of the answers to real solutions 
due to truncation errors still exist. 

To investigate the convergence question, the 480 MHz 
sweep was rerun several times with halved velocity step sizes. 
In each case the time step size was halved along with the 
velocity step size. The time step size was checked to 
insure that it was in accordance with equation 88. The 
time step size was independently varied to insure that the 
only changes seen were due to the change in the velocity 
step size. Table III summarizes the results of this invest¬ 
igation. Included are times that the computer took to run 
the case. As can be seen, the run times quickly become very 
long. Most values are not affected by changes in the velocity 
step size. The value most sensitive to the step size is the 
peak value of the velocity distribution function. Table III 
shows the percentage change from one step size to the next 
for this value. The smallest change seen was 23 $. All other 
values show percentage changes that are smaller. These small 
changes, along with long run times at small step sizes, 
convinced the author that acceptable limits of numerical 
accuracy have been reached at w* * 1 * *"/*<. and , 



































TABLE III 

Effect of Halving Velocity Increment on Numerical Results 






















































The changes in the answers are attributed to trun¬ 
cation errors. The FWHM of the force function is v* 14 m/sec. 
With a step size of 1.5 m/sec there are only about nine velocity 
steps within this FWKK. This causes truncation errors due 
to the course nature of the force function. The stability 
of the calculations insures that these errors do not 
propagate and cause the results to wildly diverge. 

However, the errors are slow to converge, and so affect 
the final results. As the velocity step size is reduced, 

the calculated answers converge to the real solutions. A 
similar argument holds for the peak of the velocity dist¬ 
ribution function. The peak becomes very narrow as the atoms 
slow down. Thus, small velocity step sizes must be used so 
enough sample points are obtained far the peak. Too large 
a velocity step size will cause the peak to not have enough 
sample points. Thus errors will be created in the calculations 
unless the velocity step size is small enough. Small enough 
will be defined as there being 20 points under the narrowest 
function. 

The final choice for the velocity step size is .375 m/sec. 
The time step size is 1.25»^ <> *" , sec. This time step size 
is in agreement with equation 88. This velocity step size 
is used since it provides ~20 points under the narrowest 
function. The narrowest function, and the function most 
sensitive to changes in the velocity step size, is the peak 
of the velocity distribution function. It changes by 2yf> 
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when the velocity step size is changed by halving it. 

If this change is neglected then all numerical results 
agree with the experimental results within 10 $. Neglecting 
the change in the peak of the velocity distribution function 
is allowable since it is only important if the peak decreses 
as the atoms are slowed. Since changes in the velocity step 
size always increase the peak value of the velocity distribution 
function, neglecting it for purposes of determining accuracy 
is acceptable. 

The important parts of this discussion are now 
summarized. The numerical results agree with the experimental 
results to within 10Ja for all values except for the FWHM 
and the peak of the distribution function. These two values 
are very sensitive to changes in velocity step size because 
of truncation error. For a step size indicated before the 
FWHM are believed to be correct to The minimum FWHM 

observed numerically agree with a calculation for the minmum 
allowable FWHM. The guassian nature of the laser is shown 
not to effect the FWHM as long as the laser is not focussed 
to micron size. The numerical calculations give rise to 
the prediction that for a given amount of frequency change, 
a faster sweep rate produces a narrower FWHM. 












































IV Conclusions and Recommendations 
Review of Ma.ior Results 

Ehrenfests theorem, along with an assumed form of the 
external electromagnetic field, yields the set of coupled 
differential equations 7 thru 10. These equations give the 
interaction force between the laser and the atoms. The in¬ 
ternal dynamics of the atom is described by the optical- 
Bloch equations 8 thru 10. This set of equations are consis¬ 
tent with a more general set of equations that are developed 
by Cook (Ref 13 pg 1087). 

In the smooth field approximation, the set of equations 
condenses to a single partial differential equation known as 
the Fokker-Planck equation. This equation is the basis for 
the analysis performed in this thesis. Along with the Fokker- 
Planck equation, coefficients are developed for the force 
and diffusion terms. These are presented in equations 15 and 
29 respectively. 

The Fokker-Planck equation is rewritten as a difference 
equation to allow numerical solutions to be obtained. Several 
difference equations are analyzed using convergence and 
stability criteria. Equation 60 is selected as the difference 
equation of choice because it has the least stringent require¬ 
ments on the step size for time and velocity. 


This difference equation is used to analyze the 
Fokker-Planck equation for a one dimensional geometry. A 
plane wave laser is incident upon a one dimensional beam. 























The laser frequency is swept at a rate given by equation 58. 
This sweeping brings different velocity atoms into resonance 
with the laser. The laser imparts quanta of momentum to the 
atoms and thus slows them down. The atomic velocity distri¬ 
bution function reflects this slowing of atoms by changing 
from a thermal distribution to a sharply peaked distribution. 

In general terms, the peak of this distribution follows the 
laser frequency and slows to low translational velocities, 
with corresponding narrow spreads in the peak. 

This analysis predicts that the translational velocities 
can be brought arbitrarily close to 0 m/sec. The minimum FWHM 
of the peak in the distribution function is 6.25 m/sec. This 
corresponds to a temperature of 1 .96 *>©''' *K, assuming the 

peak of the distribution function relaxes to a Gaussian distri¬ 
bution. This temperature is a measure of the spread in energy 
of the atoms within this slow peak. It is also a direct measure 
of the width of the peak. The minimum FWHM agrees with a limit 
developed by Cook for the smallest allowable FWHM. 

The numerical calculations suffer from truncation'error 
unless the velocity step size is small enough. Small enough 
is defined as there being 20 velocity points on the narrowest 
peak in the calculations. The narrowest peak is the distri¬ 
bution function after being cooled by the laser. The value 
of the maximum of this peak is the value most sensitive to 
changes in the velocity step size. A velocity step size of 
.1875 m/sec gives good results but the run time on the computer 

















is 8.7 CPU hours. This run time is clearly excessive. Thus, 
a velocity step size of .375 m/sec, and a time step size of 
1.25 ito"* sec are used. This results in run times of 2.5 
CPU hours. The numerical results are believed to be correct 
to within 10>1 for these step sizes. 

Conclusions 

The Fokker-Planck equation is a powerful tool to 
study the effects of resonant radiation upon atoms. All 
effects experimentally observed except one are predicted 
using a simple one dimensional analysis of the problem. The 
one dimensional imdel predicts that the atoms can be brought 
arbitrarily close to 0 m/sec in translational speed. 

Although truncation errors force a small velocity step size, 
and hence long run times, the overall method is still useful. 

The Fokker-Planck analysis predicts that for fixed frequency 
sweeps, a faster sweep rate will yield narrower FVJKK at the 
end of the interaction. 

The Fokker-Planck equation can be used for analysis 
of other types of resonant radiation interaction between 
atoms and photons. The standing wave laser used as an atomic 
trap is a gtod example. As long as the atoms enter the trap 
at slow velocities, the Fokker-Planck equation adequately 
describes their motion. 

In cases where the smooth field approximation is 
violated, the Fokker-Planck equation should be replaced by 
the quasiclassical equations developed by Cook (Ref 13 pg IO87). 
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An example of this case is found in the optical Stern-Gerlagh 
effect. Due to the quantized nature of the electric dipole 
moment, the atomic distribution function will be split in 
half by oppositely directed forces. The Fokker-Planck equation 
does not yield this result because inherent in its derivation 
is the assumption that the wave packet under consideration 
is small compared to changes in the external field. 

Recommendations 

Further work in the area of numerical analysis of 
resonate radiation can be accomplished in at least three 
different areas. No doubt other interesting applications are 
present, but these three seem the most obvious or potentially 
usefull to the author. 

The Fokker-Planck equation should be used to analyze 
the standing wave laser used as a neutral atom trap. Such 
a trap was proposed by Ashkin in 1978 (Ref 26). This type 
of trap will allow doppler free spectroscopy to be performed 
upon the atoms inside. The nature of the Fokker-Planck equation 
reverts to its general form, as in equation 112. 

|b -£ H 

The force will be given by equation 113 (Ref 13 pg 1089) 

























The diffusion coefficient is 

IH* ftU* 4 _ 


(.4 fv- -t a a. 


The Fokker-Planck equation should he analyzed in two 
dimensions for the case under consideration in this thesis. 
The form of the external field is changed to a Gaussian beam, 
as in equation 109. The force is then given by equations 49 
and 110. This analysis may allow determination of the cause 
of the observed decrease in the peak of the distribution func 
tion as the beam is slowed. 

The third area also analyzes a Guassian laser 
interacting with an atomic beam. In this case the laser is 
used to sense the position of a neutral beam. The laser 
is directed towards the neutral beam from the side. The 
atoms will be resonant with the laser for only a small angle 
of deflection from the normal. The effect of the laser 
on the emittance, and hence divergence, of the neutral 
beam would be of intrest to accelerator designers. Cooks 
quasiclassical equation would have to be solved since the 
interaction time is very short. 
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The FPLK program is used to solve the difference 
equation form of the Fokker-Planck equation. The program 
is developed to allow maximum flexibility in terms on 
analyzing different problems. This is accomplished using 
a large initial data input section. The program is based 
on a one dimensional analysis described in detail in 
Section III. The program is written in ANSI standard 
Fortran 77. 

The input data section requires that the operator 
of the program input 21 values that describe the problem 
to be analyzed. All inputed values are in the MKS system 
of units. Figure 17 shows a typical set of input values 
which were used to analyze the 480 MHz sweep discussed in 
the Analysis section. Included in figure 17 are lettered 
captions on the input values, and a key that gives a description 
of the value and its units. 

Once the input data section is mastered, almost any 
problem can be analyzed using the program. Stability 
and convergence criteria limit the step sizes for ^ and w . 

The rate of detuning must conform to equation 58. 

The CPU time required to run the program becomes 
excessive as and *»u become small. This is illustrated 
in table III, which gives run times for various cases. 

To improve execution speed, several modifications were made 
to the program. The listing given in the end of this 
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The time increments that have data printed for them 





























appendix is the modified form of the program. The changes 
amount to removing all constants which are inside the time 
loop in the program and calculating them outside the time 
loop. This saves amny calculations and speeds the program 
up by 20$. The program listed has been run for several 
different cases discussed in the thesis and yields the 
same answers. 

Figure 18 is an example of the output format of the 
program. The velocity is printed along with the value of 
the distribution function and the force at that velocity. 

This is accomplished for all velocities for one time step. 

The computer then prints the number TOTNUM (the area under 
the distribution function), and then starts the calculations 
for the next time increment. Only six values of time are 
printed in the output file since it becomes quite large. 

The file wrtdatl is used to take values from the program 
for these six times and input them into a graphics package. 
This package produces figures similar to those used in this 
thesis. The values in wrtdatl are not labled as the values 
of outputx are. The first column in wrtdatl is the velocity. 
The second column is the value of the distribution function 
at that velocity. The third is the force at that velocity 
and the fourth is the value of the time increment. 

The size of the files wrtdatl and outputx can become 
quite large for values of small velocity step size. For 
a velocity step size of .325 m/sec, the outputx file is 
























two megabytes. The wrtdatl file is JOfo smaller for this 
case. This assumes that all the values calculated are written 
into the file. The use of the variable DFS in the input 
data allows only every Ith value of the velocity to be 
printed into the outputx and wrtdatl files. Here I is the 
value assigned to the variable DFS. This can substantially 
reduce the size of the file. However, this will rule 
out being able to find the Fi/HM of the peak of the distribution 
function as it slows down. This is because the F'/JKM becomes 
very small and printing every Ith value of velocity will 
skip over the entire peak. If the F.ifKM is not desired, then 
the value for DFS can be set to allow a point every 20 m/sec 
to be printed and decent graphs will be produced. 

The graphics package used to create the figures is 

called S (Ref 27), The values in wrtdatl are read into a 

matrix in the S library. They are split into four vectors 

called force,df (value of the distribution function), vec 

(value of the velocity), and timex (value of the time increment 

the calculation is performed for). The figures generated 

for this thesis are plots of df versus vec. Thus, this is 

% 

a graph of the distribution function for time J as a function 
of velocity. Figure 19 shows the program used to plot the 
distribution function. The plotter used is a HP 7220A 
plotter. It can be used interactively with the S package. 



















































The program in figure 19 can be executed from the 
shell of the Unix system. Unix is the operating system for 
the Vax computer. The execution statement is, 

S BATCH odata pltdat 

Here odata is the name of the program in figure 19 and pltdat 
is the name of a file that results and error messages will 
appear in when the program is executed, l.'ote that S is a 
case sensitive language, and that the command S BATCH must 
be in capital letters. The ease of creating graphics programs 
in S, and the ability to execute them while in the Unix shell, 
outweighs the slow input/output times associated with the 
two software designs. The author found the combination of 
S and Unix to be a powerful one. 

The program that plots the graphics has to be executed 
while in the S library. To enter S, type S while in 
the Unix system. To exit S type q . This command 
will automatically save all files created while in S and 
return the user to the Unix shell. An on line documentation 
for 3 is available by typing help while in the S library. 

A complete description of the plotting routines is given in 
the manual for the S package (Ref 27 pg 2-16). 
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THIS PROGRAM CALCULATES THE EFFECT OF A LASER ON AN ATOMIC 
BEAM. IT IS A DIFFRENCE EQUATION OF THE FOKKER PLANCK EQUATION 
DEVELOPED FOR THIS CASE BY M G MCHARG 

REAL FORCE < 20000),DFTN< 20000),DELDFN(200000),MASS 

REAL DELTIM,DELV,TOTNUM,OMEGAC,VEL(20000),ALPHA,TEMP,A 

REAL VELO, VELMIN,F’RENUM, TIME,UK, R, VELCTR(20000> ,FPRIME< 20000) 

REAL DIFUSN(20000),DEFCTV(20000), PFTN(20000),DPRIME< 20000) 

INTEGER DFS,MT,N,Q,JCNTR 

THIS SECTION OPENS THE FILES wrtdatl AND save, wrtdatl IS THE FILE FOP 
DATA TO BE INPUTED INTO THE GRAPHICS. save IS A FILE THAT SAVES THE 
VALUE FOR JCNTR AND THE VECTORS VELCTR DFTN FORCE AND DIFUSN. 

OPEN<10,FILE =/ wrtdatl'> 

REWIND 10 

QPEN<11,FILE='save') 

REWIND 11 


PRINT*,'VARIABLE LIST' 

PRINT#,'FORCEt 2400 ELEMENT VECTOR. FORCE DISTRIBUTION FTN* 

PRINT*,'DFTN: 2400 ELEMENT VECTOR. DISTRIBUTION FTN TO BE' 

PRINT*,' USED IN BOLTZMAN EQTN.■ 

PRINT*,'MASS*. PARTICLE MASS,KILOGRAMS• 

PRINT#,'DELTIMJ INCREMENT OF TIME IN DIFFERENCING EQTN.' 

PRINT*,'DELV: CHANGE IN VELOCITY IN DIFFERENCING EQTN.' 

PRINT*,'TOTNUMJ THE AREA UNDER THE DFTN CURVE.' 

PRINT*,'VELCTR: VELOCITY ' 

PRINT*,'ALPHA: WIDTH OF DFTN IN METERS PER SECOND' 

PRINT*,'TEMPITEMPERATURE IN DEGREES KELVIN OF OVEN PRODUCING THE' 
PRINT*,' ATOMS' 

PRINT*,'A: SPONTANEOUS EMMISSION COEFFICIENT, INVERSE SECONDS' 
PRINT*,'OMEGAC: CAPITAL OMEGA, RABI FREQUENCY,INVERSE SECONDS' 
PRINT*,'VELMIN: MINIMUM VELOCITY IN DFTN.' 

PRINT#,'VELO: THE INITIAL VELOCITY WHICH FORCE AFFECTS ' 

PRINT*,'TIME: A CONSTANT USED TO TELL INTERACTION TIME' 

PRINT*,*DFS: SPACING INTEGER FOR OUTPUT PRINTING ' 

PRINT*,'PRENUM: THE PRELIMINARY TOTNUM' 

PRINT#,'MT: TOTAL NUMBER OF TIME ITERATIONS THE PROGRAM RUNS' 
PRINT*,'R.'RATE LASER FREQ IS SWEPT, IN INVERSE SECONDS' 

PRINT#,' SQUARED' 

PRINT*,'DETUNF: THE DUTUNING FREQUENCY BETWEEN LASER AND* 

PRINT#,' THE ABSORBING LINE OF THE MEDIA' 

PRINT*,'WK: WAVENUMBER. 2*PI/LAMBDA IN INVERSE METERS' 

PRINT*,'LAMBDA: THE WAVELENGTH OF THE LASER IN METERS' 

PRINT*,'DELFTN: THE CHANGE IN THE DFTN' 

PRINT#,'N: THE END NUMBER IN DO LOOPS'* 

PRINT*,'KBOLT: BOLTSMAN CONSTANT IN JOULES PER DEGREE KELVIN* 
PRINT*,'DIFUSN:20000 ELEMENT VECTOR,DIFUSION COEFFICIENT • 

PRINT*,*DEFCTV:20000 ELEMENT VECTOR,DELEFECTIVE ' 
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PRINT*,*W THE FREQUENCY OF LASER IN INVERSE SECONDS(RADIAN)* 
PRINT*,•SPDLITJ THE SPEED OF LIGHT IN NKS* 

PRINT*,'DELTIMJ TIME STEP SIZE IN SECONDS* 

PRINT*,*DELVt VELOCITY STEP SIZE IN METERS PER SECOND* 

PRINT*,’QUITJ A COUNTER SET EQUAL TO 0,1,2 FOR START,END.RESTART' 
PRINT*,* RESPECTIVELY* 

PRINT*,*JCNTR: VALUE FOR THE TIME DO LOOP. SET IT EQUAL TO 1 IF* 
PRINT*,* FIRST START,SET EQUAL TO LAST VALUE OF J IF RESTART* 
PRINT*,*B,G,D,E,FJ VALUES OF TIME COUNTER AT WHICH THE PROGRAM* 
PRINT*,* PRINTS THE OUTPUT VALUES* 

PRINT*,*HBAR: PLANCNS CONTANT DIVIDED BY 2*PI,JOULES-SECONDS* 

PRINT*,*Ct MOMENTUM OF PHOTONS IN KILOGRAM-METERS PER SECOND* 
PRINT*,'QJCOUNTER OF THE TIME INCREMENTS. WHEN IT REACHES 100* 
PRINT*,* VALUE FOR JCNTR AND THE VECTORS VELCTR,DFTN,FORCE * 
PRINT*,* AND DIFUSN ARE SAVED IN' THE SAVE FILE* 


PRINT*,* * 

THIS SECTION READS IN INITIAL VALUES 

PRINT*,* INPUT MASS,DELTIM,DELV,VELMIN,OMEGAC,TEMP,WK,MT,R,N* 

READ*, MASS,DELTIM,DELV,VELMIN,OMEGAC,TEMP,WK,MT,R,N 
PRINT*,*MASS=*,MASS,*DELTIM=*,DELTIM,*DELV=*,DELV 
PRINT*,*VELMIN=*,VELMIN,*OMEGAC=*,OMEGAC,*TEMP =*,TEMP 
PRINT*,*WK=* »WK,"MT=*»MT,*R=*,R,*N=*,N 
PRINT*,'INPUT QUIT AND JCNTR* 

READ*, QUIT,JCNTR 

PRINT*,*QUIT=*,QUIT,*JCNTR=*,JCNTR 
PRINT*,* INPUT VELO*,'INPUT DFS*,*INPUT A* 

READ*, VELO,DFS,A,W 

PRINT*,* VELO= *,VELO,* DFS* *,DFS,* A= *,A,* W= *,W 
READ*, B,G,D,E,F 

PRINT*, *B=*,B,*G=*,G,*D=*,D,*E=*,E,*F=*,F 

QUIT IS THE START/QUIT/RESTART OPTION. IF GUIT=1 THEN PROGRAM ENDS 
IF QUIT=2 THEN PROGRAM RESTARTS,IF QUIT=0 THEN PROGRAM STARTS FOR 
FIRST TIME. 

IF <QUIT .EG. 1) THEN 
GO TO 92 
END IF 

IF (QUIT .EO. 2) THEN 
READ<11,*)JCNTR 
DO 40 1=1,N 

READ <11,*)VELCTR(I),DFTN <I),FORCE(I),DIFUSN <I) 

CONTINUE 
REWIND 11 
GO TO 82 
END IF 

THIS DO LOOP INITIALIZES THE VELOCITY DISTRIBUTION FTN 
B0LTK=1.38E-23 

ALPHA=SORT< 2.*TEMP*BOLTK/MASS) 

DO 50 1=1,N 

VEL<I)=DELV*<1-1)-VELMIN 
IF<VEL<I).LT.O) THEN 
VEL<I)=0 

END IF 


V V.V V - ■ V V ■- V 







THE VELCTR SERVES TO PRINT THE REAL VELOCITY AXIS 
VELCTR<I)=DELV*< I-D-VELMIN 
BEFCTV<I)=-<WK*VELO)+WK*VELCTR<I) 

DFTN<I) = < <VEL<I)**3>/<ALPHA**4>)*EXP<-<VEL<I>**2>/<ALPHA**2>) 
CONTINUE 
PRINT*,* » 

C=WK*1.0546E-34 
HBAR=1.054E-34 
SPDLIT=2*998E8 

THE DO 70 LOOP CALCULATES THE TOTAL NUMBER TOTNUM 
PRINT*,* * 

THE FOLLOWING ARE CONSTANTS USED IN THE CALCULATIONS OF FORCE 
AND DIFUSN. 

W=SPDLIT*WK 

DIFC01=HBAR*W*0MEGAC/MASS 

DIFCO=(DIFC01/SPDLIT)*(DIFC01/SPDLIT > *< A/5») 

DIFCl=A*A+2.*OMEGAC*OMEGAC 
DIFC2=-A*A+2.*0MEGAC*0MEGAC 
DIFC3=2.*OMEGAC*OMEGAC*A*HBAR 
DIFC4=0MEGAC*0MEGAC*WK*WK*HBAR/2» 

DELDC1=DELTIM/DELV/DELV 

CINT1=<DELV-<DELTIM*R/WK>> 

DELDC2=<-BELTIM/< MASS*DELV)) 

THE DO 86 LOOP MAKES THE DIFFUSION COEFFICIENT SO THE DELDFN 
CAN BE COMPUTED. 

DO 86 1=1,N 

DIFUSN < I > = <DIFC4* < DIFC3*<12.*DEFCTV(I> +DIFC2)+A*(4.*DEFCTV(I)* 

C DEFCTV(I)+A*A)*HBAR*< 4.*DEFCTV(I)*DEFCTV(I)+A*A))/(MASS* < 4.* 

C DEFCTV <I)*DEFCTV<I)+DIFC1))/<MASS* < 4.*DEFCTV(I)*DEFCTV<I) + 

C DIFC1))/(4.*DEFCTV<I)*DEFCTV(I)+DIFC1)> +DIFCO/(4.*DEFCTV(I)* 

C DEFCTV<I)+DIFC1) 

CONTINUE 

THE DO 85 LOOP MAKES THE FORCE COEFFICIENTS SO THE DO 90 LOOP 
CAN STEP OUT IN TIME. 

DO 85 1=1,N 

FORCE <I)=-(A*OMEGAC*C*OMEGAC > / < 4 .*DEFCTV <I)*DEFCTV <I)+A#A+ 

C 2*0MEGAC*0MEGAC> 

CONTINUE 

THE DO 87 LOOP PRINTS THE INITIAL VALUES OF EVERYTHING 
PRINT*,* TIME=0 * 

DO 87 1=1,N,DFS 
PFTN<I)=DFTN(I)*100. 

PRINT*,*VEL=*,VELCTR<I),*DFTN=*,PFTN<I),*FORCE=*,FORCE(I) 
WRITE(10»990)VELCTR<I),PFTN(I),FORCE <I>,J 
FORMAT <5X,E12»5,2X,£12«5,2X,E12,5,2X,I5) 

CONTINUE 
T0TNUM=0 
DO 70 1=1,N 
PRENUM=0 

PRENUM=DELV*DFTN(I) 

TOTNUM=TOTNUM+PRENUM 

CONTINUE 

PRINT*,*TOTMUM=*,TOTNUM 







PRINT*,* * 
PRINT*,* * 


C THIS IS THE TIME LOOP. IT CALCULATES VALUES FOR THE JTH 

C ITERATION USING VALUES FROM THE J-1TH ITERATION 

82 DO 90 J=JCNTR,MT 
TIME=J*DELTIM 

C THE DO 91 LOOP INTERPOLATES TO MOVE THE FORCE AND DIFFUSAION 

C COEFFICIENTS DOWN THE VELOCITY AXIS. THIS KEEPS THE PROGRAM 

C FROM RECALCULATING THE COEFFICIENTS EACH TIME. 

DO 91 1=1,N-l 

DPRIME < I ) =-(CINT1* < DIFUSN(I+1>-DIFUSN <I))/DELV)+DIFUSN(1 + 1) 
FPRIMECI)=-<CINT1*<FORCE<1+1)-FORCE(I)>/DELV>+FORCE(1+1) 

91 CONTINUE 

C THE DO 100 LOOP FINISHES MOVING THE FORCE AND DIFFUSION 

C COEFFICIENTS DOWN THE VELOCITY AXIS. 

DO 100 1=1,N 
DIFUSN<I)=DPRIME<I) 

FORCE <I)=FPRIME <I) 

100 CONTINUE 

C THE DO 80 LOOP CHANGES THE DISTRIBUTION FTN FOR VELOCITY 

DO 80 I=2,N-l 

DELDFN<I)=DELDC2*(FORCE(I+1)*DFTN<I+1) 

C -FORCE(I)*DFTN<I))+<DELDCl*<DIFUSN(1+1>*DFTN(1+1)-2»* 

C DIFUSN<I)*DFTN<I)+DIFUSN<1-1)*DFTN(1-1))> 

80 CONTINUE 

C THE DO 81 LOOP CHANGES THE DISTRIBUTION FUNCTION FOR VELOCITY 

DO 81 1=1,N 

DFTN(I)=DFTN(I)+DELDFN(I) 

81 CONTINUE 
• 

C THE DO 75 LOOP PRINTS THE VELOCITY,DFTN,FORCE. 

IF<J.EQ.B.OR.J.EQ.G.OR.J.EQ.D.OR.J.EQ.E.OR.J.EQ.F)THEN 
PRINT*,*TIME=*,TIME 
DO 75 1=1,N,DFS 

PFTN<I)=BFTN<I)*100 

PRINT*,‘VELOCITY**,VELCTR<I),*DFTN=*,PFTN<I>,*FORCE=*,FORCE(I) 
WRITEUO, 1000) VELCTR< I ) ,PFTN< I ) ,FORCE < I), J 
1000 FORMAT(5X,E12.5,2X,E12»5,2X,E12.5,2X,I5) 

75 CONTINUE 

PRINT#,*TOTNUM=*,TOTNUM 
END IF 


C THE DO 76 LOOP RECALCULATES TOTNUM 

TQTNUM=0 
DO 76 L=1,N 
PRENUM=0 

PRENUM=DELV*DFTN< L) 

TOTNUM®TOTNUM+PRENUM 
76 CONTINUE 






JCNTR=J 

THE Q COUNTER PRINTS EVERY 100 TIME AND TOTNUM 
Q=Q+1 

IFCG.EG,100)THEN 
Q=0 

WRITE<11,*)JCNTR 
DO 77 1=1,N 

WRITE<11,#)VELCTR<I),DFTN<I),FORCE(I) »DIFUSN<I) 

CONTINUE 

REWIND 11 

PRINT*,*TIME=',TIME 
PRINT*,' ' 

PRINT*,''T0TNIJM=' ,TOTNUM 

PRINT*,' ' 

PRINT*,' * 

PRINT*,' ' 

END IF 
CONTINUE 
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